Otherware/Otherware_1_SB_Development.iso/mac/misc/biology/treedraw.sit/Tree Draw Deck/card_9880.txt

Text File | 1990-06-26 | 26.0 KB | 969 lines

-- card: 9880 from stack: in
-- bmap block id: 0
-- flags: 0000
-- background id: 9703
-- name: PlotGramXCMD.p


-- part contents for background part 3
----- text -----
DrawGramXCMD.p

-- part contents for background part 2
----- text -----
(* drawgramXCMD.p
 *
 * Mac Hypercard XCMD
 * Usage:
 *     drawgram tree, [parameters,] [output field]
 *       returns HP drawing in The Result
 *
 *       tree == container w/ tree description
 *       param. == cont. w/ parameters
 *       output fld == bkgnd field for text messages, if desired
 *
 * hacked by d.g.gilbert from DrawGram.pas of Phylip package
 * dogStar Software && Indiana University Biology Dept.
 * email: [e-mail]
 *
 * Language: MPW-Pascal 3.0
 *

 #-------- Makefile
  pascal drawgramXCMD.p -mbg off
  Link -w -rt XCMD=10041 Γêé
    -m ENTRYPOINT Γêé
    -sg drawgram Γêé
    drawgramXCMD.p.o Γêé
    "{Libraries}HyperXLib.o" Γêé
    "{Libraries}Interface.o" Γêé
    "{PLibraries}PasLib.o" Γêé
    "{PLibraries}SANElib.o" Γêé
    "{Libraries}Runtime.o" Γêé
    -o "Tree Draw Deck"
 #------------------------
 *
 *
 *)


{$R-}
{$Z+} (* has to have this or else put EntryPoint in Interface section *)

UNIT DummyUnit;


INTERFACE

USES
  Types,
  SANE,
  QuickDraw,
  Memory,
  OSUtils,
  ToolUtils,
  Packages,
  PasLibIntf,
  HyperXCmd;


IMPLEMENTATION


(* XCMD entry -- MUST BE FIRST CODE IN SEG *)

procedure drawgram(outfile, treefile, parmfile, plotfile: ptr);
  forward;

PROCEDURE EntryPoint( pXCmd: XCmdPtr);
CONST
    maxOut  = 2000;
    maxPlot = 70000; {?? this should be expandable, but we'll be lazy}
    Times   = 20; {Font unit}
VAR
    hTree, hParams, hOut, hPlot: Handle;
    outfld  : str255;
    wind    : grafPtr; sfont, ssize : integer;

BEGIN

  hTree := pXCmd^.params[1];
  if hTree = NIL then begin
    pXCmd^.returnValue := pasToZero(pXcmd,
            '§§§ Tree description is missing ');
    exit( EntryPoint);
    end;
  HLock(hTree);

  hParams:= pXCmd^.params[2];
  if (hParams = NIL) then hParams:= NewHandleClear(10);
  HLock(hParams);

  hOut  := pXCmd^.params[3];
  if hOut = NIL then outfld:= ''
  else zerotopas( pXCmd, hout^, outfld);
  hOut  := NewHandleClear(maxOut);
  HLock(hOut);

  hPlot := NewHandleClear(maxPlot);
  HLock(hPlot);

  getPort(wind);
  sfont:= wind^.txFont; ssize:= wind^.txSize;
  TextFont( Times); TextSize( 12); {set for label size calcs}

  DrawGram(hOut^, hTree^, hParams^, hPlot^);

  HUnlock(hPlot);
  HUnlock(hOut);
  HUnlock(hParams);
  HUnlock(hTree);
  TextFont( sfont); TextSize( ssize);

  if outfld <> '' then SetFieldByName( pXCmd, false, outfld, hOut);
  DisposHandle( hOut);

  pXCmd^.returnValue:= hPlot;
END;


function min(a,b:longint):longint;
begin
  if (a<b) then min:= a else min:= b;
end;

function max(a,b:longint):longint;
begin
  if (a>b) then max:= a else max:= b;
end;


FUNCTION Strtod(s: str255; var ends: integer): real;
(* stdlib strtod, strtol use %GlobalData !!  *)
var
  x   : real;
  drec: decimal;
  valid: boolean;
BEGIN
  {- ends := 1; -- let caller set this.  1 is lowest for pstr}
  if ends < 1 then ends:= 1;
  str2dec( s, ends, drec, valid);
  StrToD := dec2num( drec);
END;


{ hypercard won't allow globals in XCmds that are needed
  by standard file handlers }

type  carray = packed array [0..32000] of char;
      cp  = ^carray;

function hcEof( f: ptr): boolean;
begin
  hcEof:= cp(f)^[0] = chr(0);
end; {hcEof}

function  hcEoln( f: ptr): boolean;
begin
  hcEoln:= (cp(f)^[0] = chr(0))
        or (cp(f)^[0] = chr(13));
end; {hcEoln}


procedure hcRead( var f: ptr; var c: char);
begin
  c:= cp(f)^[0];
  if c <> chr(0) then longint(f):= longint(f) + 1;
end; {hcRead}

procedure hcReadReal( var f: ptr; var x: real);
var
  i, e : integer;
  s   : str255;
begin
 i:= 0;
 while (cp(f)^[i] >= ' ') and (i < 255) do begin
   s[i+1]:= cp(f)^[i];
   i:= i+1;
   end;
 s[0]:= chr(i);
 e := 1;
 x := Strtod( s, e);
 longint(f):= longint(f) + e;
end; {hcReadReal}

procedure hcReadInt( var f: ptr; var i: integer);
var  x: real;
begin
  hcReadReal(f, x);
  i:= trunc(x);
end; {hcReadInt}

procedure hcReadLn( var f: ptr);
begin
  while (cp(f)^[0]<> chr(13))
    and (cp(f)^[0]<> chr( 0)) do
      longint(f):= longint(f) + 1;
  if (cp(f)^[0] = chr(13)) then
      longint(f):= longint(f) + 1;
end; {hcReadLn}

procedure hcWrite( var f: ptr; s: str255);
var  j, l: longint;
begin
  l:= length(s);
  for j:= 1 to l do cp(f)^[j-1] := s[j];
  cp(f)^[l+1]:= chr(0);
  longint(f):= longint(f) + l;
end; {hcWrite}

procedure hcWriteInt( var f: ptr; n: integer);
var  s: str255;
begin
  numtostring( n, s);
  hcWrite( f, s);
end;

procedure hcWriteReal( var f: ptr; x: real; w,d: integer);
var  s  : decstr;
     decfm: decform;
begin
   decfm.style := fixedDecimal;
   decfm.digits:= d;
   Num2Str(decfm, x, s);
   while length(s) < w do s:= concat(' ',s);
   hcWrite( f, str255(s));
end; {hcWriteReal}

procedure hcWriteln( var f: ptr; s: str255);
var  s1: string[1];
     outh : handle;
     lout : longint;
begin
  s1 := '1'; s1[1]:= chr(13);
  hcWrite(f, concat(s, s1));
end; {hcWriteln}


function QDfontHeight: integer;
var  fi: fontInfo;
begin
  getFontInfo(fi);
  QDfontHeight:= fi.ascent + fi.descent + fi.leading;
end;





procedure drawgram(outfile, treefile, parmfile, plotfile: ptr);
  (* Version 3.3.  Copyright (c) 1986, 1990 by Joseph Felsenstein and
    Christopher A. Meacham.
    Permission is granted to copy, distribute, and modify this
    program provided that (1) this copyright message is not removed
    and (2) no fee is charged for this program.

    hacked w/ a macheté by d.gilbert, June 1990:
    -- dropped a plotters but hp
    -- replaced NEW with NewPtr and DisposeTree
    -- dropped all font stuff
    -- replaced all file i/o with pointer i/o
    -- dropped all interactive stuff
       and replaced w/ parmfile reader
    -- parmfile format:
      -- 1 space b/n param name_value, uppercase is minimum abbrev.
      Grow Horizontal/Vertical
      Style Cladogram/Phenogram/Eurogram/Swoopogram/curVogram
      Uselengths Yes/No
      Rotate 95.0
      Depth 1.23
      Length 1.23
      Nodelength 1.23
      Position Intermediate/Weighted/Centered/iNner/Vshaped

    *)

  CONST
    maxnodes = 300;
    maxnch = 30;
    point = '.';
    pi = 3.141592653;
    epsilon = 0.00001;

    kGrows  = 'G';      {parameter keys, for parmfile}
    kStyle  = 'S';
    kUseLengths = 'U';
    kLabelRot = 'R';
    kTreeDepth = 'D';
    kStemLength = 'L';
    kNodeSpace = 'N';
    kNodePosition = 'P';

    HersheyfontHeight = 28;

  TYPE
    growth = (vertical, horizontal);
    treestyle = (cladogram, phenogram, curvogram, eurogram, swoopogram);
    penstatustype = (penup, pendown);
    plotstring = packed ARRAY[1..maxnch] OF CHAR;
    nodeptr = ^node;
    node = RECORD
      next, back : nodeptr;
      tip : BOOLEAN;
      nayme : plotstring;
      naymlength, tipsabove : INTEGER;
      xcoord, ycoord, level : REAL
      END;

  VAR
    ntips, nextnode, nextnext,
    numtochange, oldx, oldy, nmoves, payge : INTEGER;
    haslengths, uselengths : BOOLEAN;
    xmargin, ymargin, topoflabels, rightoflabels, leftoflabels, tipspacing,
      scale, xscale, yscale, xoffset, yoffset, nodespace, stemlength,
      treedepth, xnow, ynow, xunitspercm, yunitspercm, xsize, ysize, xcorner,
      ycorner, oldxhigh, oldyhigh, oldylow, oldxlow, treeline, labelline,
      labelheight, labelrotation,  expand, rooty : REAL;
    grows : growth;
    style : treestyle;
    root : nodeptr;
    nodep, nextp : ARRAY[1..maxnodes] OF nodeptr;
    penchange, oldpenchange : (yes, no);
    nodeposition : (weighted, intermediate, centered, inner, vshaped);


  PROCEDURE uppercase (VAR ch : CHAR);
  BEGIN
    IF ((ch >= 'a')  AND (ch <= 'z'))
    THEN ch := CHR(ORD(ch) + ORD('A') - ORD('a'));
  END; (* uppercase *)

  PROCEDURE treeread;
    (* read a tree from the treefile and set up nodes and pointers *)
    VAR ch : CHAR;

    PROCEDURE getch (VAR c : CHAR);
      (* get next nonblank character *)
    BEGIN
      REPEAT
        IF hcEOLn(treefile) THEN hcReadLN(treefile);
        hcRead(treefile, c);
      UNTIL c <> ' ';
    END; (* getch *)

    PROCEDURE addelement (VAR p : nodeptr; q : nodeptr);
      (* read in and add next part of tree, it will be node p
        and will be hooked to pointer q *)
      VAR pfirst : nodeptr;
        n : INTEGER;
        notlast : BOOLEAN;

      PROCEDURE processlength(p : nodeptr);
        VAR digit, ordzero : INTEGER;
          valyew, divisor : REAL;
          pointread : BOOLEAN;
      BEGIN (* processlength *)
        ordzero := ORD('0');
        pointread := FALSE;
        valyew := 0.0;
        divisor := 1.0;
        getch(ch);
        digit := ORD(ch)-ordzero;
        WHILE ((digit >= 0) AND (digit <= 9)) OR (ch=point) DO BEGIN
          IF ch = point
          THEN pointread := TRUE
          ELSE BEGIN
            valyew := valyew*10.0 + digit;
            IF pointread
            THEN divisor := divisor*10.0;
            END;
          getch(ch);
          digit := ORD(ch)-ordzero;
          END;
        p^.level := valyew/divisor;
      END; (* processlength *)

    BEGIN (* addelement *)
      nextnode := nextnode + 1;
      ptr(p) := NewPtr(sizeof(node)); {-NEW(p);}
      nodep[nextnode] := p;
      IF ch = '(' THEN BEGIN
        p^.tip := FALSE;
        p^.tipsabove := 0;
        pfirst := p;
        notlast := TRUE;
        WHILE notlast DO BEGIN
          ptr(p^.next) := NewPtr(sizeof(node)); {- NEW(p^.next);}
          p := p^.next;
          nextnext:= nextnext + 1;
          nextp[nextnext]:= p;    {!save for dispose}
          p^.tip := FALSE;
          getch (ch);
          addelement (p^.back, p);
          pfirst^.tipsabove := pfirst^.tipsabove + p^.back^.tipsabove;
          IF ch = ')' THEN BEGIN
            notlast := FALSE;
            REPEAT  getch (ch);
            UNTIL (ch = ':') OR (ch = ',') OR (ch = ')') OR (ch = ';');
            END;
          END;
        p^.next := pfirst;
        p := pfirst;
        END
      ELSE BEGIN
        p^.tip := TRUE;
        p^.tipsabove := 1;
        ntips := ntips + 1;
        n := 1;
        REPEAT
          IF (ch = '_') THEN ch := ' ';
          IF n < maxnch
          THEN p^.nayme[n] := ch;
          IF hcEOLn(treefile) THEN hcReadLN(treefile);
          hcRead(treefile, ch);
          n := n + 1;
        UNTIL ((ch = ':') OR (ch = ',') OR (ch = ')'));
        IF n > maxnch
        THEN n := maxnch + 1;
        p^.naymlength := n - 1;
        END;
      IF ch = ':'
      THEN processlength (p)
      ELSE haslengths := haslengths AND (q = NIL);
      p^.back := q;
    END; (* addelement *)

  BEGIN (* treeread *)
    haslengths := TRUE;
    ntips := 0;
    nextnode := 0; nextnext:= 0;
    getch (ch);
    addelement (root, NIL);
    hcReadLN(treefile);
    uselengths := haslengths;
  END; (* treeread *)


  procedure disposeTree;
  { release Ptrs for Mac/Hypercard }
  var  i: integer;
  begin
    for i:= 1 to nextnode do DisposPtr(ptr(nodep[i]));
    for i:= 1 to nextnext do DisposPtr(ptr(nextp[i]));
  end;


  PROCEDURE plotrparms;
    (* set up initial characteristics of plotter or printer *)
  BEGIN
    xcorner := 0.0;
    ycorner := 0.0;
{Hewlett-Packard plot setup}
    penchange := yes;
    xunitspercm := 400.0;
    yunitspercm := 400.0;
    xsize := 24.0;
    ysize := 18.0;
  END; (* plotrparms *)

  PROCEDURE initialparms;
  BEGIN
    plotrparms;
    xmargin := 0.08 * xsize;
    ymargin := 0.08 * ysize;
    xscale := xunitspercm;
    yscale := yunitspercm;
    style := cladogram;
    grows := vertical;
    labelrotation := 45.0;
    nodespace := 3.0;
    stemlength := 0.05;
    treedepth := 0.5/0.95;
    IF uselengths
    THEN nodeposition := intermediate
    ELSE nodeposition := vshaped;
  END; (* initialparms *)

  PROCEDURE getparms;
    (* get from user the relevant parameters for the plotter and diagram *)
  VAR key, ch : CHAR;
      ok : BOOLEAN;
      x : real;
  BEGIN
  while not hcEOF(parmfile) do begin
    {read parameter key/name}
    repeat hcRead(parmfile, key); until (key > ' ') or hcEOLn(parmfile);
    {skip rest of parameter name...}
    repeat hcRead(parmfile, ch);  until (ch <= ' ') or hcEOLn(parmfile);

    uppercase(key);
    CASE key OF

      kGrows: BEGIN  {! note change in documentation }
        hcread(parmfile, ch);
        uppercase(ch);
        case ch of
          'H' : grows := horizontal;
          'V' : grows := vertical;
          end;
        END;

      kStyle: BEGIN
        hcread(parmfile, ch);
        uppercase (ch);
        CASE ch OF
          'C' : style := cladogram;
          'P' : style := phenogram;
          'E' : style := eurogram;
          'S' : style := swoopogram;
          'V' : style := curvogram
          END;
        END;

      kUseLengths: BEGIN
        hcread(parmfile, ch);
        uppercase (ch);
        if haslengths then CASE ch OF
          'Y': begin uselengths:= true;
               nodeposition:= intermediate;
               end;
          'N': begin uselengths:= false;
               nodeposition:= vshaped;
               end;
            end;
        END;

      kLabelRot: BEGIN
        hcREADReal(parmfile, labelrotation);
        END;

      kTreeDepth: BEGIN
        hcReadReal(parmfile, treedepth);
        END;

      kStemLength: BEGIN
        hcReadReal(parmfile, x);
        if (x >= 0.0) AND (x < 0.9) then stemlength := x;
        END;

      kNodeSpace: BEGIN
        hcReadReal(parmfile, x);
        if (x <> 0.0) then nodespace := 1.0/x;
        END;

      kNodePosition: BEGIN
        hcRead(parmfile, ch);
        uppercase(ch);
        CASE ch OF
          'W' : nodeposition := weighted;
          'I' : nodeposition := intermediate;
          'C' : nodeposition := centered;
          'N' : nodeposition := inner;
          'V' : nodeposition := vshaped
          END;
        END;

      END;
    hcReadln(parmfile);
    end;
  END; (* getparms *)



  PROCEDURE calculate;
    (* compute coordinates for tree *)
    VAR sum, tipx, maxtextlength, textlength,
      firstlet, maxheight, fontheight, angle : REAL;
      i : INTEGER;

    PROCEDURE calctraverse (p : nodeptr; lengthsum : REAL);
      (* traverse to establish initial node coordinates *)
      VAR x1, y1, x2, y2, x3, w1, w2, sumwx, sumw, nodeheight : REAL;
        pp, plast : nodeptr;
    BEGIN (* calctraverse *)
      IF p = root
      THEN nodeheight := 0.0
      ELSE IF uselengths
        THEN nodeheight := lengthsum + p^.level
        ELSE nodeheight := 1.0;
      IF nodeheight > maxheight THEN maxheight := nodeheight;
      IF p^.tip
      THEN BEGIN
        p^.xcoord := tipx;
        IF uselengths
        THEN p^.ycoord := nodeheight
        ELSE p^.ycoord := 1.0;
        tipx := tipx + tipspacing;
        END
      ELSE BEGIN
        sumwx := 0.0;
        sumw := 0.0;
        pp := p^.next;
        x3 := 0.0;
        REPEAT
          calctraverse (pp^.back, nodeheight);
          sumw := sumw + pp^.back^.tipsabove;
          sumwx := sumwx + pp^.back^.tipsabove*pp^.back^.xcoord;
          IF ABS(pp^.back^.xcoord-0.5) < ABS(x3-0.5)
          THEN x3 := pp^.back^.xcoord;
          plast := pp;
          pp := pp^.next;
        UNTIL pp = p;
        x1 := p^.next^.back^.xcoord;
        x2 := plast^.back^.xcoord;
        y1 := p^.next^.back^.ycoord;
        y2 := plast^.back^.ycoord;
        CASE nodeposition OF
          weighted : BEGIN
            w1 := y1 - nodeheight;
            w2 := y2 - nodeheight;
            IF (w1 + w2) <= 0.0
            THEN p^.xcoord := (x1 + x2)/2.0
            ELSE p^.xcoord := (w2*x1 + w1*x2)/(w1+w2);
            END;
          intermediate : p^.xcoord := (x1 + x2)/2.0;
          centered : p^.xcoord := sumwx/sumw;
          inner : p^.xcoord := x3;
          vshaped : p^.xcoord := (x1 + x2 + (y1 - y2)/maxheight)/2.0
          END;
        IF uselengths
        THEN p^.ycoord := nodeheight
        ELSE BEGIN
          p^.ycoord := (x1 - x2 + y1 + y2)/2.0;
          IF nodeposition = inner
          THEN BEGIN
            IF ABS(x1-0.5) > ABS(x2 - 0.5)
            THEN BEGIN
              p^.ycoord := y1 + x1 - x2;
              w1 := y2 - p^.ycoord;
              END
            ELSE BEGIN
              p^.ycoord := y2 + x1 - x2;
              w1 := y1 - p^.ycoord;
              END;
            IF w1 < epsilon
            THEN p^.ycoord := p^.ycoord - ABS(x1-x2);
            END;
          END;
        END;
    END; (* traverse *)

    FUNCTION lengthtext(pstring : plotstring; nchars : INTEGER) : REAL;
       VAR i, j, code : INTEGER;
          cfix, sumlength, heightfont, widthfont : REAL;
    BEGIN
       sumlength := 0.0;
       heightfont:= HersheyfontHeight;
       cfix:= HersheyfontHeight / QDfontHeight;
       FOR i := 1 TO nchars DO BEGIN
          widthfont:= round( cfix * charwidth(pstring[i]));
          sumlength := sumlength + widthfont;
          END;
       lengthtext := sumlength;
    END; (* lengthtext *)

  BEGIN (* calculate *)
    maxheight := 0.0;
    maxtextlength := 0.0;
    IF nodep[1]^.naymlength > 0
    THEN firstlet := lengthtext (nodep[1]^.nayme, 1)
    ELSE firstlet := 0.0;
    sum := 0.0;
    tipx := 0.0;
    FOR i := 1 TO nextnode DO
      IF nodep[i]^.tip
      THEN BEGIN
        textlength := lengthtext (nodep[i]^.nayme, nodep[i]^.naymlength);
        IF textlength > maxtextlength
        THEN maxtextlength := textlength;
        END;
    fontheight := HersheyfontHeight;
    angle := pi*labelrotation/180.0;
    maxtextlength := maxtextlength/fontheight;
    textlength := textlength/fontheight;
    firstlet := firstlet/fontheight;
    IF ntips > 1
    THEN labelheight := 1.0/(nodespace*(ntips-1))
    ELSE labelheight := 1.0/nodespace;
    IF angle < (pi/6.0)
    THEN tipspacing := (nodespace
                  + COS(angle)*(maxtextlength-0.5))*labelheight
    ELSE IF ntips > 1
      THEN tipspacing := 1.0/(ntips-1.0)
      ELSE tipspacing := 1.0;
    topoflabels := labelheight*(1.0 + SIN(angle)*(maxtextlength-0.5)
                          + COS(angle)*0.5);
    rightoflabels := labelheight*(COS(angle)*(textlength-0.5) +SIN(angle)*0.5);
    leftoflabels := labelheight*(COS(angle)*firstlet*0.5+SIN(angle)*0.5);
    calctraverse (root, sum);
    rooty := root^.ycoord;
    FOR i := 1 TO nextnode DO
      nodep[i]^.ycoord := stemlength*treedepth + (1.0-stemlength)*treedepth
         *(nodep[i]^.ycoord-rooty)/(maxheight-rooty);
    rooty := 0.0;
  END; (* calculate *)

  PROCEDURE rescale;
    (* compute coordinates of tree for plot or preview device *)
    VAR i : INTEGER;
      treeheight, treewidth, extrax, extray, temp : REAL;
  BEGIN (* rescale *)
    treeheight := 0.0;
    FOR i := 1 TO nextnode DO
      IF nodep[i]^.ycoord > treeheight
      THEN treeheight := nodep[i]^.ycoord;
    treeheight := treeheight + topoflabels;
    treewidth := (ntips-1)*tipspacing + rightoflabels + leftoflabels;
    IF grows = vertical
    THEN BEGIN
      expand := (xsize - 2*xmargin)/treewidth;
      IF (ysize - 2*ymargin)/treeheight < expand
      THEN expand := (ysize - 2*ymargin)/treeheight;
      extrax := (xsize - 2*xmargin - treewidth*expand)/2.0;
      extray := (ysize - 2*ymargin - treeheight*expand)/2.0;
      END
    ELSE BEGIN
      expand := (ysize - 2*ymargin)/treewidth;
      IF (xsize - 2*xmargin)/treeheight < expand
      THEN expand := (xsize - 2*xmargin)/treeheight;
      extrax := (xsize - 2*xmargin - treeheight*expand)/2.0;
      extray := (ysize - 2*ymargin - treewidth*expand)/2.0;
      END;
    FOR i := 1 TO nextnode DO BEGIN
      nodep[i]^.xcoord := expand*(nodep[i]^.xcoord + leftoflabels);
      nodep[i]^.ycoord := expand*(nodep[i]^.ycoord - rooty);
      IF grows = horizontal
      THEN BEGIN
        temp := nodep[i]^.ycoord;
        nodep[i]^.ycoord := expand*treewidth-nodep[i]^.xcoord;
        nodep[i]^.xcoord := temp;
        END;
      nodep[i]^.xcoord := nodep[i]^.xcoord + xmargin + extrax;
      nodep[i]^.ycoord := nodep[i]^.ycoord + ymargin + extray;
      END;
    IF grows = vertical
    THEN rooty := ymargin + extray
    ELSE rooty := xmargin + extrax;
  END; (* rescale *)

  PROCEDURE plot(pen : penstatustype; xabs, yabs : REAL); {for HP}
  BEGIN
    IF pen=pendown THEN hcWrite(plotfile, 'PD')
                   ELSE hcWrite(plotfile, 'PU');
    hcWriteInt(plotfile, ROUND(xabs));
    hcWrite(plotfile, ',');
    hcWriteInt(plotfile, ROUND(yabs));
    hcWriteLN(plotfile, ';');
  END; (* plot *)

  PROCEDURE plottree (p, q : nodeptr);
    (* plot part or all of tree on the plotting device *)
    CONST segments = 40;
    VAR i : INTEGER;
      x1, y1, x2, y2, x3, y3, f, g, h, fract, minny, miny : REAL;
      pp : nodeptr;
  BEGIN (* plottree *)
    x2 := xscale*(xoffset + p^.xcoord);
    y2 := yscale*(yoffset + p^.ycoord);
    IF p <> root
    THEN BEGIN
      x1 := xscale*(xoffset + q^.xcoord);
      y1 := yscale*(yoffset + q^.ycoord);
      plot (penup, x1, y1);
      CASE style OF
        cladogram : plot (pendown, x2, y2);
        phenogram : BEGIN
          IF grows = vertical
          THEN plot (pendown, x2, y1)
          ELSE plot (pendown, x1, y2);
          plot (pendown, x2, y2);
          END;
        curvogram : FOR i := 1 TO segments DO BEGIN
            f := i/segments;
            g := i/segments;
            h := 1.0 - SQRT(1.0-g*g);
            IF grows = vertical
            THEN BEGIN
              x3 := x1*(1.0-f)+x2*f;
              y3 := y1 + (y2 - y1)*h;
              END
            ELSE BEGIN
              x3 := x1 + (x2 - x1)*h;
              y3 := y1*(1.0-f)+y2*f;
              END;
            plot (pendown, x3, y3);
            END;
        eurogram : BEGIN
          IF grows = vertical
          THEN plot (pendown, x2, (2*y1+y2)/3)
          ELSE plot (pendown, (2*x1+x2)/3, y2);
          plot (pendown, x2, y2);
          END;
        swoopogram : IF (grows = vertical) AND (NOT (ABS(y1-y2) < epsilon))
              OR ((grows = horizontal) AND (NOT (ABS(x1-x2) < epsilon)))
          THEN BEGIN
            IF grows = vertical
            THEN miny := p^.ycoord
            ELSE miny := p^.xcoord;
            pp := q^.next;
            WHILE pp <> q DO BEGIN
              IF grows = vertical
              THEN minny := pp^.back^.ycoord
              ELSE minny := pp^.back^.xcoord;
              IF minny < miny
              THEN miny := minny;
              pp := pp^.next;
              END;
            IF grows = vertical
            THEN miny := yscale*(yoffset+miny)
            ELSE miny := xscale*(xoffset+miny);
            IF grows = vertical
            THEN fract := 0.3333*(miny-y1)/(y2-y1)
            ELSE fract := 0.3333*(miny-x1)/(x2-x1);
            FOR i := 1 TO segments DO BEGIN
              f := i/segments;
              IF f < fract
              THEN g := f/fract
              ELSE g := (f-fract)/(1.0-fract);
              IF f < fract
              THEN h := fract*SQRT(1.0-(1.0-g)*(1.0-g))
              ELSE h := fract + (1.0-fract)*(1.000001 - SQRT(1.000001-g*g));
              IF grows = vertical
              THEN BEGIN
                x3 := x1*(1.0-f)+x2*f;
                y3 := y1 + (y2 - y1)*h;
                END
              ELSE BEGIN
                x3 := x1 + (x2 - x1)*h;
                y3 := y1*(1.0-f)+y2*f;
                END;
              plot (pendown, x3, y3);
              END;
            END
        END;
      END
    ELSE BEGIN
      IF grows = vertical
      THEN BEGIN
        x1 := xscale*(xoffset + p^.xcoord);
        y1 := yscale*(yoffset + rooty);
        END
      ELSE BEGIN
        x1 := xscale*(xoffset + rooty);
        y1 := yscale*(yoffset + p^.ycoord);
        END;
      plot (penup, x1, y1);
      plot (pendown, x2, y2);
      END;
    IF NOT p^.tip
    THEN BEGIN
      pp := p^.next;
      WHILE pp <> p DO BEGIN
        plottree (pp^.back, p);
        pp := pp^.next;
        END;
      END;
  END; (* plottree *)

  PROCEDURE plotlabels;
    TYPE pentype = (treepen, labelpen);
    VAR i : INTEGER;
      compr, dx, dy, angle : REAL;
      lp : nodeptr;

    PROCEDURE changepen(pen : pentype);
    BEGIN (* changepen *)
      CASE pen OF
        treepen:  hcWriteLN(plotfile, 'SP1;');
        labelpen: hcWriteLN(plotfile, 'SP2;');
        END;
    END; (* changepen *)

  PROCEDURE plottext(pstring : plotstring; nchars : INTEGER;
    height, compress,  x, y, slope : REAL);
  CONST xstart = 10;
     ystart = 35;
  VAR i, j : INTEGER;
     sinslope, cosslope : REAL;
     penstatus : penstatustype;
     s1 : string[1];
  BEGIN (* plottext *)
    sinslope := SIN(pi*slope/180.0);
    cosslope := COS(pi*slope/180.0);
    plot( penup, x, y);
    hcWrite(plotfile,'DI ');
    hcWriteReal(plotfile, cosslope,7,4);
    hcWrite(plotfile,', ');
    hcWriteReal(plotfile, sinslope,7,4);
    hcWriteln(plotfile,';');
    hcWrite(plotfile,'LB');
    for i:= 1 to nchars do hcWrite(plotfile,pstring[i]);
    s1:= '1'; s1[1]:= chr(3);
    hcWriteln(plotfile,s1);
  END; (* plottext *)

  BEGIN (* plotlabels *)
    compr := xunitspercm/yunitspercm;
    IF penchange = yes THEN changepen(labelpen);
    angle := labelrotation*pi/180.0;
    FOR i := 1 TO nextnode DO BEGIN
      IF nodep[i]^.tip
      THEN BEGIN
        lp := nodep[i];
        dx := -labelheight*expand*0.70710*COS(angle+pi/4.0);
        dy := labelheight*expand*(1.0-0.70710*SIN(angle+pi/4.0));
        IF grows = vertical
        THEN plottext(lp^.nayme,lp^.naymlength,
          labelheight*expand*xscale/compr, compr,
          xscale*(lp^.xcoord+dx+xoffset), yscale*(lp^.ycoord+dy+yoffset),
          -labelrotation)
        ELSE plottext(lp^.nayme,lp^.naymlength,
          labelheight*expand*yscale, compr, xscale*(lp^.xcoord+dy+xoffset),
          yscale*(lp^.ycoord-dx+yoffset), -labelrotation+90.0)
        END;
      END;
    IF penchange = yes THEN changepen(treepen);
  END; (* plotlabels *)

  PROCEDURE initplotter;
  BEGIN
    hcWriteLN(plotfile, 'IN;SP1;VS10.0;');
  END;

  PROCEDURE finishplotter;
  BEGIN
    plot(penup, 1.0, 1.0);
    hcWriteLN(plotfile,'SP;AF;');
  END;

  PROCEDURE drawit;
  BEGIN
    xoffset := 0.0;
    yoffset := 0.0;
    plottree (root, root);
    plotlabels;
  END; (* drawit *)


BEGIN (* drawgram *)
  hcWriteLN(outfile,'DRAWGRAM from PHYLIP version 3.3');
  hcWriteLN(outfile,'Reading tree ... ');
  treeread;
  hcWritelN(outfile,'Tree has been read.');
  initialparms;
  getparms;
  calculate;
  rescale;
  initplotter;
  hcWriteLN(outfile,'Writing plot file ...');
  drawit;
  finishplotter;
  disposeTree;
  hcWriteLN(outfile,'Finished.');
END; {drawGram}

END. {DummyUnit}